function lpcvd_1DTRISnapshot_plot(filename,LatticeSize,Thickness,Mode)
fid = fopen(filename);
data = textscan(fid,'%f %f','HeaderLines',3);
fclose(fid);
layer = data{1}; % NOTE: the value in layer starts from 0
site = data{2};  % NOTE: the value in site is from 0 to (LatticeSize-1)

TopLayer = max(layer);
if nargin == 2
    BottomLayer = 0;
    Thickness   = TopLayer;
else
    BottomLayer = TopLayer-Thickness+1;
end
indx = layer>=BottomLayer;
layer = layer(indx);
site = site(indx);
%% Use RECTANGLE command

key = [((BottomLayer)*LatticeSize):(LatticeSize*max(layer)+LatticeSize-1)]';
key_vac = (layer)*LatticeSize+site;
j = 1;
flag_occ = ones(size(key));
num_vac = length(key_vac);
for i = 1:length(key)
    if key(i)>key_vac(j)
        error('Something wrong');
    elseif key(i)==key_vac(j)
        flag_occ(i) = 0;
        j = j+1;
    end
    if j > num_vac
        break;
    end
end
key_occ = key(flag_occ == 1);
x = mod(key_occ,LatticeSize);
y = (key_occ-x)./LatticeSize;
for i = 1:length(key_occ)
   if mod(y(i),2)==1
        x(i) = x(i)+0.5;
   end
end
box on;
for i = 1:length(key_occ)
    rectangle('position',[x(i)-0.5,y(i)+0.15,1,1],'FaceColor',[0.90,0.90,0.90],'EdgeColor',[0.50,0.50,0.50],'Curvature',[1,1]);
end
YRange = ylim();
ylim([min(y),min(y)+Thickness*1.2]);
xlim([-0.5,LatticeSize]);

% scatter(x,y);
% ylim([BottomLayer,TopLayer+0.2*(TopLayer-BottomLayer)]);

